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ABSTRACT 

Wc present an algorithm for performing precise aperture photometry on criticaUy 
sampled astrophysical images. The method is intended to overcome the small-aperture 
limitations imposed by point-sampling. Aperture fluxes are numerically integrated 
over the desired aperture, with sinc-interpolation used to reconstruct values between 
pixel centers. Direct integration over the aperture is computationally intensive, but 
the integrals in question are shown to be convolution integrals and can be computed 
> 10000 X faster as products in the wave- number domain. The method works equally 
well for annular and elliptical apertures and could be adapted for any geometry. A 
sample of code is provided to demonstrate the method. 



1 INTRODUCTION 



In aperture photometry, one seeks to measure the flux which 
would have been received had the light passed through a 
physical aperture (e.g. a photomultipher tube). An aper- 
ture flux in an astrophysical image is typically measured as 
the sum of the pixel fluxes for pixels near a given source. 
The region summed over is the 'aperture', and may contain 
whole and/or partial or weighted pixels. Circular apertures 
are common for astrophysical measurements, and are imple- 
mented in some form in mos t standard soft ware packages , 
including IRAF dTod vl Il986l), DAOPHO T (IStetson|[l987l ). 
and SExtractor (Ber lin fc Arnouts|[l996l ). 

Aperture fluxes are most commonly used in calibration 
when constructing a 'growth curve'. A growth curve plots 
the enclosed flux within an aperture as a function of its ra- 
dius, and is essential for calibrating fluxes measured through 
an optimal (typically small) aperture to bright isolated stan- 
dard stars which must be measured with a large aperture 
to include all flux. Also, a number of commonly used flux 
measurements are aperture fluxes computed for scientiflcally 
interes ting ape rture shapes or sizes. Examples include Kron 
fluxes l|Kron|[l980l ) and Petrosian fluxes ( Petrosian 1976 ) . 

For aperture radii comparable to the pixel dimensions, 
if no partial pixels are used, the 'circular' aperture can be 
quite non-circular (e.g. alxlor2x2 square, a cross, 
etc.). As the flux may change signiflcantly across a pixel, 
the use of geometric partial pixels (e.g. including 0.5 x flux 
for a pixel with half of its area enclosed in the aperture) 
will not accurately represent the flux contained in the corre- 
sponding fraction of the pixel. To maintain a truly constant 
aperture geometry (circular or otherwise), the flux must be 
interpolated between pixels and integrated. 

Here we describe a straightforward aperture photom- 
etry algorithm, which under normal observing conditions 
preserves circular aperture geometry perfectly and remains 
well-defined to arbitrarily small radii. For the method to be 



accurate, the pixels must critically sample the point-spread 
function (PSF); i.e. the spatial frequency of the pixels must 
be at least twice the highest spatial frequency present in 
the PSF. The method works equally well for annular and 
elliptical apertures, and could be adapted for any desired 
shape. The ability to handle elliptical apertures, in par- 
ticular, makes it useful for computing Kron and elliptical 
Petrosian fluxes. It is implemented in the image processing 
pipeline under dev elopment for the Large Synoptic Survey 
Telescope (LSST) l|lvezic et al.ll2008l) . 

To compute the flux that would have been received 
through a truly circular aperture centred on a point source, 
we perform a continuous reconstruction of the discretely 
sampled pixel fluxes (Section [2)) and integrate to yield the 
circular aperture flux (Section[5)|. The two-dimensional inte- 
gral is computationally intensive, but we use Fourier meth- 
ods to compute it quickly without loss of precision (Sec- 
tion 13]). In Section \5\ we modify the method to compute 
annular and elliptical apertures; and in Section [6] we evalu- 
ate the performance of the algorithm by comparing growth 
curves for model point-spread functions (PSFs) to their 
known values. Results are discussed and summarized in Sec- 
tions [7] and |S1 respectively. Some additional calculations are 
included in Appendix [Xj and an example program (written 
in Python) is included in Appendix IB] 



2 RECONSTRUCTING A DISCRETELY 
SAMPLED PSF 

A telescope produces a continuous function and delivers it 
to a detector, which convolves with the pixel response and 
point-samples it. In our technique, we use the sampling theo- 
rem to recover the continuous function, and measure it. Our 
only assumption is that the function must be band-limited. 
In this context, this implies that the Fourier transform of the 
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Figure 1. The sine interpolation of a one-dimensional PSF. A 
discretely sampled PSF (a) is transformed with a discrete Fourier 
transform (b), and truncated at the Nyquist wave number (c). 
This truncation is the product of a boxcar and the wave-number 
coefficients, and is therefore represented by the convolution of a 
sine (the Fourier transform of a boxcar) with the sampled points 
in real space (d). If the sampled function is truly band limited, 
the sine representation shown as a dashed line in (d) is exact. 



PSF has negligible po-wer above the Nyquist wave number 
which corresponds to the pixel sampling. 

The sampling theorem states that a discretely sampled 
continuous function can be completely represented by its 
samples. In one dimension 



sin(7r(a; — Xi)) 

■k{x — Xi) 



(1) 



Equation [T] can be used to interpolate p{x) given the 
discrete samples pi. As the sin(a;)/a; term is a sine function, 
this is known as sine interpolation. Figure [T] demonstrates 
how the sine interpolation of a one-dimensional PSF can 
be understood in terms of Fourier methods. If the sampled 
function is truly band-limited (no power above the Nyquist 
wave number), the sine interpolation is an exact represen- 
tation of the underlying continuous function. The proof of 
the samp l ing th eorem can be found in the original work by 
IShannonI l| 19491 ). 

The sine decomposition of a two-dimensional function 
has the form 



p{x,y) ^^Yl 



Pij 



sin(7r(a; - a;^)) sin(7r(j/ - j/^)) 



n{x - 



(2) 



2.1 Testing the Band-Limit of a Gaussian PSF 

The requirements for the sine reconstruction to be valid are 
that the PSF be band limited, and that an infinite sequence 
of values are available. However, in both cases more practical 
limits are sufficient. Here we test the first assumption (the 
band limit) in the context of reconstructing a continuous 
PSF. Requirements on the extent of the data are dealt with 
in Section r3.il 

The validity of the band-limit assumption was verified 
by evaluating the error in integrated flux as a function of 
the sampling interval. This was estimated analytically with 
a circular bivariate single- Gaussian PSF model, and also 
measured numerically with a double-Gaussian PSF. A PSF 
is typically well-modeled by a double Gaussian (a sum of 
two bivariate Gaussians, one modeling the PSF core and 
the other modeling the wings). The single-Gaussian model 
used for the analytic estimate is intended to test the core 
Gaussian of a double- Gaussian PSF. If the 'band limit' cri- 
terion is satisfied for the narrower Gaussian, it would easily 
be satisfied for the larger one representing the wings. 

For our analytic model, the calculation of the frac- 
tional error due to aliased power in a sampled Gaussian is 
a straightforward exercise, and is presented in Appendix 1X1 
The resulting upper-limit estimate for the fractional flux er- 
ror is 

error _ e"'' " 
total flux ^ 27ri/2cr ' ^ ' 

To verify the method more rigorously, it was also tested 
numerically with known PSFs. For the numerical tests, the 
so-called 'double Gaussian' is a reasonable approximation to 
a typical observed PSF. Our double-Gaussian test PSF has 
core and wing components, with the wing component having 
10% of the core's amplitude and a width of o"wing = 2(Tcorc- 
A double Gaussian also has the advantage that it can be 
integrated analytically to provide an exact known flux for 
comparison to the sine-integrated value, and it was therefore 
used for these numerical tests. 

The response of the pixel is considered to be a part of 
the PSF structure (i.e. the continuous PSF produced by the 
telescope is convolved with a pixel and then point-sampled 
by the detector) , and the double Gaussian model used in the 
numerical tests was directly point-sampled. As convolving 
a Gaussian with a pixel is very close to convolving with a 
another Gaussian having cr = 1/ \/T^, the resulting function 
remains very close to a Gaussian. 

Both theoretical (equation|3]) and numerically measured 
errors are shown as a function of the PSF size in pixels (i.e., 
the sampling) in Figure [2] An aperture with a radius of 6 
pixels was used. For the numerical tests, sub-pixel centroid 
shifts were applied to verify the method for PSF profiles that 
are not centred at integer pixel positions, and representative 
curves corresponding to shifts of {5x,5y) = (0.0,0.0); and 
(0.5, 0.5) are presented. Provided the PSF's a^orc is larger 
than ~1 pixel (corresponding to a FWHM > 2.5 pixels for 
a double Gaussian), fractional flux errors due to undersam- 
pling are <^ 0.001, or ~1 mmag. 

It is tempting to assume that any error introduced 



The RMS width of a unit boxcar function (i.e. a pixel) is l/y/vl. 
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Figure 2. The flux error due to truncation/aliasing of wave num- 
bers with (J > TT as a function of PSF size in pixels for a double 
Gaussian (lower panel), with the absolute value shown in the up- 
per panel. The flux error is 1— measured/true, where 'measured' 
refers to our sinc-integration algorithm (aperture radius of 6 pix- 
els), and 'true' refers to the exact integral of the double Gaussian 
(sum of circular bivariate Gaussians having different amplitudes 
(wing being 10% of core) and widths (awing = 2(Tcorc))- Curves 
are shown for a double Gaussian centred on an integer pixel value 
(solid), and one shifted by 0.5 pixels (dashed). Limits based on 
our conservative theoretical estimate are shown with dot-dashed 
lines. The theoretical estimate is the limit for truncated flux, while 
the sinc-integrated curves include aliased flux. Dotted lines indi- 
cate a fractional error of 0.001, or ~1 mmag. For the Gaussian 
PSFs, the flux error is <C Immag provided sampling is c > 1 
pixel. In cases of severe undersampling, the formalism of includ- 
ing the pixel response in the double-Gaussian PSF breaks down; 
therefore, values for small PSF widths (ccoro < 0.4 pixels) are 
not shown. 

through ahasing would be systematic and correspond to a 
fractional excess or deficit of fiux. Such an error could be cor- 
rected by calibrating the measured fiuxes against standard 
stars. ISowever, Figure [2] shows that the errors are differ- 
ent for the {S^,5y) = (0.0,0.0); and (0.5,0.5) pixel shifts. In 
general, the error for each measured source depends on the 
sub-pixel shift of the PSF and cannot be corrected through 
calibration. Further, galaxies and other non-point-sources 
will be composed of different spatial frequencies and their 
errors will be different from those of point sources and from 
other galaxies. 

If the method is to be used in the construction of a 
growth curve (flux measurements at a series of progressively 
larger apertures), the band-limited nature of the PSF must 
be taken seriously as the error is a function of radius and 
can be large for undersampled PSFs. For additional informa- 
tion, see Section [6.11 where the performance of the method 
is tested directly by measuring growth curves. 

With these caveats in mind, a double-Gaussian PSF is 
band limited for the pixel scales and aperture sizes com- 
monly used in astronomy (cr > 1 pixel, or FWHM > 2.5 
pixels) . 



3 COMPUTING APERTURE FLUX 

With the discrete pixel fluxes, Pij, interpolated to form a 
continuous function, p(x,y) (see equation [2]) , a precise aper- 
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Figure 3. The Wij coefficients shown as images (bottom), and 
cross-sections through the middle row (top) for apertures with 
radii r = 1 pix (left), r = 3 pix (middle), and r = 10 pix (right). 



ture flux, /a, can be obtained by integrating over the desired 
aperture A: 

/a = J J p{x,y) dxdy. (4) 

A 

The integration can be performed to arbitrary preci- 
sion, but depending on the structure of the aperture, and 
integrator used, it may be computationally intensive (recall 
that the p{x, y) function is a sum of sine functions - one for 
each pixel). 

Fortunately, the Pij terms in equation [5] are constant 
and can be factored out of the integral in equation [l] allow- 
ing the aperture flux to be expressed as a weighted sum: 

j i 

where the weighting terms are 

A 

Thus, for each pixel we compute a corresponding 
weight, Wij, in the sum for the aperture flux. The integral 
(wij) does not depend on pij. It can be computed once and 
shifted (see Section 14. ip to be applied to different sources 
in a given frame. This is critical for efficiency. Without the 
ability to pre-compute the Wij coefficients, several numerical 
integrals would need to be computed for each source being 
measured, and the method would be prohibitively slow. 

For a large circular aperture, pixels near the center have 
weight Wij ~ 1. Near the edge of the aperture, values transi- 
tion from ~1 to ~0. However, values outside the aperture do 
not go completely to 0. Some non-zero weight remains even 
several pixels outside the aperture radius. Examples of Wij 
coefficients for different-sized circular apertures are shown 
as images and cross-sections in Figure [S] 
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Footprint width: N^Jr^ 



Figure 4. The flux error occuring at different footprint sizes 
(lower panel). Absolute values are shown on a semilog scale (up- 
per panel) to allow small values to be distinguished. The flux error 
for a well-sampled (cr = 1.2 pixel) double-Gaussian PSF is influ- 
enced by the width of the footprint used to compute the aperture 
flux. The pixels outside the aperture carry weight, and the foot- 
print must be larger than the aperture to include them. Results 
for three aperture radii are shown: r=1.0, 2.0, and 5.0. Footprint 
widths are normalized by the aperture radius. Footprint widths 
smaller than two aperture radii would explicitly exclude a por- 
tion of flux and were not considered. The error is large only when 
small apertures are used with small footprints. A footprint hav- 
ing a 4 pixel border around the aperture is sufficient to achieve 
< 0.001 fractional flux error for even the smallest aperture. 




Figure 5. The noise effective area (J^^ ui^/vrr^, top panel), and 
the signal to noise ratio (bottom panel) for circular sine apertures 
as functions of aperture radius in pixels. For large apertures, the 
noise effective area approaches 1 and S/N approaches that of a 
unity-weighted 'tophat' aperture (shown dashed). For small radii, 
the noise effective area is smaller for the sine aperture, and the 
S/N approaches a constant rather than zero. The S/N curve is 
based on a double-Gaussian PSF with (Tcorc=1.5 pixel. Optimal 
S/N for both methods is near ~ 1.6crcorc. 



3.1 The Minimum Size Footprint 

As the Wij coefficients are non-zero outside the aperture 
radius, they cannot be ignored and the footprint for the 
aperture (the x riy image containing the Wij coefficients) 
must be larger than the aperture itself. To determine the re- 
quired footprint size, Gaussian PSFs were planted in images 
of varying size: 2 — 12 x rap. These were compared to the 
known values for a double Gaussian analytically integrated 
over the same aperture. A well-sampled double Gaussian 
(ccorc ~ 1.2 pixels) was used with aperture radii of rap=1.0, 
2.0, and 5.0 pixels. For small apertures (< 2 pixel), errors 
can be as large as a few percent if the footprint is too small. 
For large apertures, errors are negligible even for footprints 
only slightly larger than the aperture radius. Empirically, we 
find that with a 4 pixel border around the aperture, errors 
are reduced to < 0.001 for a well-sampled PSF. Figure [4] 
shows the flux error as a function of the width of the foot- 
print. 

3.2 The Effects of Noise 

Poisson noise is always present in an astrophysical image and 
has a flat power spectrum which is not band limited. The 
sine aperture flux / can be regarded as a simple weighted 
sum of pixel intensities, li. 

i 

The variance of the sum, assuming a pixel variance af 

is 



l^j^'^WiGi. (8) 

i 

Equations [7] and |8] remain true regardless of the choice 
of weight values, and the variances of two candidate weight- 
ing methods can be compared as a ratio. For constant pixel 
noise ai (i.e. sky limited observations), the sine aperture 
compared to a top hat aperture having Wi = 1 gives a ratio 
of 

-LJ-U...-J2^Vn. (9) 

i 

The number of pixels A'^ is the area of the comparison 
apert ure, and can be regarded as a noise equivalent 

area l|Kinelll9"83 ). i.e. the area of a unity- weighted aperture 
which would contribute the same variance to the measure- 
ment. For the sine aperture, Figure [5] (upper panel) shows 
equation [5] versus aperture radius, with the denominator re- 
placed by the area of the circular aperture being measured. 

For large aperture radii, the noise equivalent area 
asymptotically approaches 1, but it decreases for small aper- 
tures. The noise in a sine fiux is less than would be ex- 
pected in a unity- weighted 'top hat' aperture! The inclusion 
of pixels outside the aperture radius might be expected to 
introduce additional noise, but instead serves to reduce the 
overall variance. The pixels are present outside the aperture 
radius, but they carry less weight and the noise is averaged 
over a larger sample of pixels. 

Note that we have not considered the shape of the PSF. 
We have demonstrated only that a measurement performed 
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with a circular sine aperture will have lower variance than 
one performed with a 'top hat' aperture. By considering the 
PSF, we can determine the signal to noise ratio as a function 
of aperture radius. This is shown for a double Gaussian with 
o"corc = 1.5 pixels in the lower panel of Figure (5] Curves for 
both sine and top hat apertures are shown. The peak S/N 
at ~ 1.6 X (Tcore ~ 2.4 pixels is a known result for a Gaussian 
PSF (see for example IPritchet fc Klin3ll98ll '). However, as 
the S/N deteriorates to zero for the top hat aperture, it 
levels off to a constant for the sine aperture. Thus, though 
S/N is not optimal at small radii, it approaches a constant 
rather than zero. 



4 FAST COMPUTATION OF COEFFICIENTS 

The Wij double integral in equation [6] can be pre-computed, 
but is still computationally intensive. However, by appeal- 
ing to Fourier methods a fast computation algorithm is 
available. By expressing the aperture as a top hat function 
with radius p, n((a;^ + y'^)^^^ /2p), and noting that the sine 
functions are even (i.e., sinc(7r(a; — Xi)) = sinc(7r(a;i — a;))), 
the integral can be expressed as the convolution of a two- 
dimensional Cartesian sine with a circular top hat, sampled 
discretely at Xi, yj-. 

J J Tv{x^-x) -^(yj-y) 

(10) 

The convolution can be computed quickly as a product 
of the Fourier-transformed functions in the wave-number do- 
main. Both functions have analytic transforms: 



Benchmark calculations of the Wij coefficients for aper- 
tures with radii r = 1 — 10 pixels found ~ 10000 x speed 
increases (< 1 ms versus several seconds) when calculations 
were performed in the wave-number domain. 



4.1 Shifting the Aperture Centroid 

The centroid of a PSF is typically measured with sub-pixel 
accuracy, and it is essential that it be possible to recenter 
our aperture on non-integer pixel coordinates. If the Wij 
coefficients are pre-computed for integer pixel coordinates, 
sinc-interpolation can be used to shift them by the necessary 
sub- pixel offset, S:,:,Sy. 

As previously discussed, sinc-interpolation is valid only 
for functions that are band limited. We demonstrated that a 
PSF can be treated as a band-limited function, but have not 
done so for the Wij components. However, the Wij compo- 
nents represent the aperture convolved with a sine function. 
As a sine function is itself band limited, it follows by the 
convolution theorem, that any convolution with it will also 
form a band-limited function. The Wij coefficients can safely 
be shifted by sinc-interpolation. 

It is also possible to apply a sub-pixel shift directly to 
Wij coefficients when they are created with the FFT method. 
To perform aperture photometry on several sources, shifting 
pre-computed Wij coefficients is faster than directly com- 
puting them for each source individually, but only slightly. 
If this approach is taken, shifting can be achieved easily in 
the wave-number domain by applying the Shift Theorem 
(/(a; -S)^ e~"'""=F(fc)). Dropping the redundant n{k/n) 
terms, equation 1131 then becomes 



n 



(^^ + y^)^ 

2p 



p Jl (2TVp{ki + kt,)2 



(11) 



and 



sine {7t{x, - x)) sine {7v{y, -y))^U (^^ j H (^^ j . (12) 

The former yields an Airy function, and the latter yields 
a two-dimensional boxcar that extends to the Nyquist wave 
numbers. The coefficients can be obtained by computing the 
product of the wave-number domain components from equa- 
tions [11] and [121 c'lid taking the inverse Fourier transform: 



w — T 



p Jl {2-Kp{kl + kl)^^^ 

(fc| -f fc^) 2 



(14) 

Due to the shift, the wave- number domain components 
are neither real- valued nor even. However, the application of 
the shift term, e~'^'^^*, preserves conjugate symmetry (i.e., 
if f(x) is real- valued, then F(k) — F*{~k), where '*' de- 
notes complex conjugation). With conjugate symmetry, the 
spatial domain values (i.e., in pixel space) are guaranteed 
to be real-valued. This is not surprising given the physical 
interpretation of a shift in the spatial domain. 



n I — n 

TT J \ TT 



ky\ P Jl (2^p(fcl + fcg)i) 



(13) 

This can be computed with a fast Fourier transform 
(FFT). The array of values in the wave- number domain ex- 
tends out to the Nyquist wave number - precisely where the 
boxcar truncates the Airy function. Thus, instead of per- 
forming the real-space integral in equation |S] the Wij coeffi- 
cients can be found directly by computing the values for the 
Airy function (equation [TTJ , and taking the inverse FFT. 
The transformed Wij coefficients are entirely real-valued, 
and the coefficients are mathematically identical to those 
shown in Figure [S] 



5 OTHER APERTURES 

Photometry with any aperture can be performed with this 
method; the circular aperture is conventional in observa- 
tional astronomy, and conveniently has an analytic form in 
the wave-number domain. If the Fourier transform of the 
desired aperture can be accurately computed, it can be sub- 
stituted for the Airy function in equations 1131 and 1141 

If the desired aperture has no analytic transform, the 
sine convolution integral in equation [6] can be computed di- 
rectly. As this is more computationally intensive, it could be 
impractical to recompute for each source being measured, 
and the integral should be pre-computed and shifted to the 
pixel coordinates of the object being measured. 



© 2013 RAS, MNRAS OOO.mfTTI 



6 S.J. Bickerton 



-^L^ JLr- tJ Ltltl 



--'nmillll |lllr'"'Ml| ||llniln~- 



■ 2 4 6 S 101214 2 4 6 8 101214 10 20 30 40 50 60 




O O 



2 4 6 8 101214 



10 20 30 40 50 60 

Pixels 



)=HI)°, e=0.h 



) = HI)°, e=l).5 



Lm 



a = ll)-lSpLx 
9 = HI)° , e=l).5 



2 4 6 8 101214 10 20 30 40 50 60 




2 4 6 8 101214 




10 20 30 40 50 60 



Figure 6. The Wij coefficients shown as images (bottom), and 
cross-sections through the middle row (top) for annular apertures 
with inner/outer radii r = 1 — 2 pix (left), r = 3 — 5 pix (middle), 
and r = 10 — f 5 pix (right). 

5.1 Annular Apertures 

Due to the addition theorem {f{x) + g{x) <4> F{k) + G{k)), 
an annular aperture can be constructed trivially by taking 
a difference of Airy terms representing the inner and outer 
radii of the annulus. The Wij coefficients for a range of an- 
nular apertures are shown in Figure [6l 



Figure 7. The Wij coefficients shown as images (bottom), and 
cross-sections through the middle row (top) for annular elliptical 
apertures with inner/outer semi-major axes a = 1 — 2 pix (left), 
a = 3 — 5 pix (middle), and a = 10 — 15 pix (right). Ellipses are 
rotated hy 6 = 30° and have ellipticity e = 1 — b/a = 0.5, where 
b and a are the semi-minor and semi-major axes, respectively. 

this work, the popular Moffat PSF (|Moffatl [19691 ') was also 
tested. In the second test, our aperture magnitudes were 
compare d to values prod uced by the SExtractor software 
package (|Bertin fc Arnouts.lQQS 'l. 



5.2 Elliptical Apertures 

The special case of an elliptical aperture can be handled 
by applying the Similarity Theorem {f{ax) <^ \a\^^F{k/a)) 
to the circular case. The compression of a circle along a 
given axis forms an ellipse, and the corresponding Fourier 
transform can be formed by dilating that of a circle along 
the appropriate axis. 

The Fourier transform preserves rotation; if f{x, y) is 
rotated by 6 degrees, F{kx,ky) will be rotated by a corre- 
sponding amount in the same sense. We therefore compress 
the Airy function along the ky axis, and rotate it by the 
desired amount with a rotation matrix applied to the kx, ky 
wave numbers. 

To demonstrate all the properties presented, Wij coeffi- 
cients for a selection of annular elliptical apertures are shown 
in Figure [T] 

In Appendix|Bl we include a sample of Python code that 
demonstrates the method by computing growth curves for 
model PSFs. The coefficient images in Figures [3l [6l and \7\ 
were constructed with the wijCoeff icientsO routine in- 
cluded in the sample code. 



6 PERFORMANCE AND EFFECTIVENESS OF 
THE ALGORITHM 

The accuracy of the method was tested by constructing and 
measuring images of sources with known PSFs and fluxes. 
Two different tests were performed. In the flrst test, aperture 
magnitudes were measured for progressively larger apertures 
to determine the growth curves for model PSFs, and these 
were compared to those of the known, analytic solutions. 
In addition to the double Gaussian PSF used throughout 



6.1 Measuring the Grovirth Curves for Double 
Gaussian and Moffat PSFs 

For simplicity, analytic double Gaussian and Moffat PSFs 
were used for testing. Neither is a band-limited function, 
and different widths ((Tcoro values for the double Gaussians 
and FWHM for the Moffat PSFs) were tested to evaluate the 
behaviour near the critical sampling limit, Ucorc = 1 pixel. 
The pixel values for the double Gaussian were computed at 
coordinates x, y with 



Pi] 



2n (a|,,, + bol- 



+ h e 







(15) 



Here, (Twing ~ 2(Tcorc, and h = 0.1. Values for the Moffat PSF 
were computed with 



13-1 



1.0 + 



2 I 2 

X +y 



(16) 



where l3 = 4.765 (|Truiillo et al.l 120011 ). and a = 

FWHM/2(2^'''^ - 1)^/^ The FWHM used to compute the 
Moffat PSF was 2.4670crcorc. Note that this differs slightly 
from the value of FWHM=2.3548ct for a single Gaussian. 

Note that these are a sampling of the PSFs as we do 
not integrate over the pixels. The PSF is considered to be a 
combination of multiple contributing components, including 
the aperture response (an Airy function), the atmospheric 
blurring, and the pixel response. 

Three widths were evaluated: (TcorG=0.8, 1.0, and 1.2 
pixels, corresponding to under-, critical-0, and well-sampled 



Strictly speaking, there is no 'critical' sampling limit for either 
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Aperture Radius [Pixels] Aperture Radius [Pixeis] 

Figure 8. Measured growth curves and residuals with respect 
to computed theoretical fluxes. Left and right panels show values 
for double-Gaussian and Moffat PSFs, respectively. PSF widths 
of (7=0.8, 1.0, and 1.2 pixels were used to evaluate under-, critical- 
, and over-sampling. The differences between the measured and 
theoretical curves are too small to be seen directly, but are shown 
as residuals in the three lower panels. Residuals shown with a 
dashed line are the minimum observed, and correspond to no 
sub-pixel shift. Those shown with a solid line are the maxima, 
and correspond to a half-pixel shift in both x and y coordinates. 
The curves are very similar in structure, but differ dramatically 
in amplitude. The method is clearly unsuitable for sub-sampled 
(cr < 1.0 pixel) PSFs, but performs extremely well for cr > 1.2 
pixel. 



PSFs. Sources were planted in the center of a 16x16 pixel 
image. Additional sub-pixel ofTsets were applied to test the 
efTects of shifting the aperture. 

Sinc-integrated aperture fluxes were measured from the 
images with aperture radii, < r < 4 pixels (0.1 pixel incre- 
ments). The known PSFs (eguation llSp were then integrated 
directly over the same apertures to obtain the true theoret- 
ical fluxes. The structure of the residuals varies depending 
on the sub-pixel shift applied. For both PSFs, the smallest 
residual was observed when no sub-pixel offset was applied, 
and the largest residual was observed when half-pixel offsets 
were applied to both x and y. Discrepancies are highest for 
apertures with radii < 1 pixel. Growth curves of the mea- 
sured fluxes are presented with residuals in Figure [S] To 
illustrate the effects of the sub-pixel shift, the curves shown 
are those representing the extrema. 

The amplitudes of the residual curves decrease rapidly 
as sampling improves. Discrepancies as high as a few percent 
can be seen for the undersampled PSFs, and the method is 
clearly unsuitable under such conditions. The well-sampled 
PSF (o-core = 1.2, FWHM= 3.0) performs extremely weU 
at all aperture sizes (fractional residual less than 0.001 and 
0.002 for double Gaussian and Moffat, respectively). 



a Gaussian or a Moffat PSF as they are not band-limited func- 
tions. We use the term loosely, but for all practical intents and 
purposes (as we show) (TcorG=1.2 pixels is quite reasonable. 



6.2 Comparison to SExtractor Aperture 
Magnitudes 

We also compared sine magnitudes to t he values produced 
by t he SExtractor software package (|Bertin fc Arnouti 
Il996l ) version 2.3.2. SExtractor is one of the most widely 
used photometry packages in the astronomy community and 
is the most suitable algorithm for such a comparison. 

Double Gaussian PSFs were added to the center of 
24x24 pixel images. Tests were performed with well-sampled 
(o"core=l-5 pixel) and under-sampled (crcore=0.7 pixel) dou- 
ble Gaussians. As with previous tests, the wing component 
was 2 X wider and had an amplitude which was 10% of the 
core component's amplitude. Each test was run with and 
without noise for aperture radii of r=1.0, 2.0, and 5.0 pix- 
els. The tests performed with noise included a 100 count 
sky level, and the pixel values were replaced with random 
Poisson variates. For simplicity, a gain of 1.0 was used. 

For a given test, 500 trials were run for fluxes cover- 
ing several magnitudes. In each case, the double Gaussian 
was given a random sub-pixel offset. The test image was 
written to a FITS file and was measured with both SEx- 
tractor and the sine algorithm. In order to ensure that only 
the photometric algorithms were compared, the sine algo- 
rithm was run using the x, y pixel centroids as measured by 
SExtractor. The PHOT_APERTURES parameter was set 
to 2x the radius (SExtractor specifies aperture diameter), 
and MAG_ZEROPOINT was set to an arbitrary value of 
25. For tests including noise and sky, the known sky level of 
100 counts was subtracted before processing, and SExtrac- 
tor's BACK.TYPE parameter was set to MANUAL to dis- 
able internal background estimation. Results of the testing 
are shown in Figures |5] (well-sampled PSF), and [TU] (under- 
sampled PSF). 

For well-sampled data measured with a large (r >5 
pixel) aperture, the performance of the two algorithms is 
similar. For smaller apertures, the sine algorithm outper- 
forms SExtractor. In the case of under-sampled data with a 
large aperture, SExtractor outperforms. However, this per- 
formance deteriorates for smaller apertures while the sine 
algorithm's does not. For r <1.0 pixel apertures, the sine 
method yields the better result. 



7 DISCUSSION 

The method we have described relies on very well- 
established mathematical theorems, namely the sampling 
theorem, and various theorems associated with the Fourier 
transform. However, various aspects of our analysis require 
some discussion. 

Our decision to use pure sine interpola tion in place of a 
Lanczos kernel (a tapered sine function, see lLanczod (|l956l )) 
requires some justification. When we integrate the sine func- 
tion over the aperture, it is defined within the entire region 
of integration, regardless of the location of the sine peak (i.e. 
the pixel in Wij being considered). The value obtained for 
each Wij pixel is therefore more accurate than that which 
would be obtained by a Lanczos, which tapers to zero be- 
yond its region of support (typically 4 to 7 pixels). How- 
ever, this difference is truly negligible. Our principal rea- 
son for using a pure sine was that it offers a convenient 
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Figure 9. Magnitude error versus known magnitude for sine 
(dots) and SExtraetor (+ symbols) aperture magnitudes of well- 
sampled PSFs. A lower sub-plot shows the standard deviation of 
points in bins with width 0.5 mag. Well-sampled double-Gaussian 
PSFs (0-00^=1.5 pixel) were added to 24 X 24 pixel images and 
were measured with both algorithms. Aperture radii of r=1.0 
(top), 2.0 (middle), and 5.0 (lower) pixels were used. Left panels 
show results for which no noise was added to the double-Gaussian 
PSFs, and right panels show results which included Poisson noise 
and the addition of a 100 count sky level. For apertures with 
r >5.0 pixels, performance is similar for both algorithms. For 
smaller apertures, the sine method outperforms in both precision 
and accuracy. 



Fourier transform (unity throughout our region of interest 
in the wave-number domain). In order to minimize demand 
on computing resources for e.g. a Kron flux, computation in 
the wave-number domain is essential. The Fourier transform 
of a Lanczos kernel is sufficiently complicated that the pure 
sine is the better choice. 

As there are a variety of PSF models available, an al- 
ternative to measuring aperture fluxes directly would be to 
fit an appropriate model and integrate the model over the 
desired aperture. When the object being measured is not a 
point source, this approach is often preferred. In such cases, 
the objects have different sizes and shapes, and use of a con- 
stant aperture radius would not allow objects to be mean- 
ingfully compared. Model fluxes can frequently be the best 
choice in such situations. However, for the purpose of cali- 
bration, the use of a model introduces an unnecessary source 
of potential systematic error. A poorly chosen model will 



Figure 10. Magnitude error versus known magnitude for sine 
(dots) and SExtraetor (+ symbols) aperture magnitudes of under- 
sampled PSFs. The panels are the same as those described in 
Figure|9] In this case, under-sampled (0.7 pixel) double- Gaussian 
PSFs were used. For large apertures the SExtraetor magnitudes 
are more precise and more accurate than the sine magnitudes. For 
small apertures, the sine magnitudes are more precise and more 
accurate, becoming the better option for apertures with r <1.0 
pixel. 

never fit properly and will corrupt the aperture correction 
and thus the zero point of the calibration. 

As a more practical verification of the method, the sinc- 
interpolated aperture fluxes are being used to compute aper- 
ture fluxes, and to determine the aperture correction for the 
LSST photometric pipeline that is currently under devel- 
opment. The sine method has been found to perform ex- 
tremely well, providing (mpsF — map) < 0.001 (verifying 
aperture correction), and (map — meat) ^ 0.001 (verifying 
measurement accuracy). Here, mpsF and map are PSF and 
(sine) aperture magnitudes, and meat is a catalog magni- 
tude corresponding to the known input flux in a simulated 
image. Several hundred simulated images have been gener- 
ated, processed, and tested, covering a variety of observ- 
ing conditions (cloud, seeing, etc). The LSST simulator uses 
a ray-tracing algorithm to trace photons through a simu- 
lated atmosphere and through simula ted optics to produce 
each test image l|Peterson et al.|[2013l . in preparation). The 
ray-tracing simulations include realistic fleld-dependent and 
wavelength-dependent effects that result in complex point- 
spread functions and photometric variation that challenge 
photometry algorithms. . The aperture magnitudes show no 
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Figure 11. Si ne aperture photomet ry performed on a simulated 
LSST i' image jPeterson et al.ll2013l . in preparation). For clarity, 
only 1/4 of the values from a full frame are shown. The upper 
panel shows map — meat versus meat , and the lower panel shows 
the standard deviation of the upper panel data in 0.5 mag bins. 
The aperture radius was 7 pixels (1.4"), and magnitudes are com- 
pared to the exact known values (i.e. meat) used as input for the 
simulator. 

systematic trends with magnitude, and the variabihty in 
shape and PSF width observed across a focal plane of in- 
dividual CCD images is handled well by a 2D polynomial 
aperture correction. An example of sine aperture photom- 
etry performed on LSST simulated data is shown in Fig- 
ure [11] As yet, no failure modes have been observed for the 
algorithm. 

As stated earlier, any error introduced as a result of 
a non-band-limited PSF will not be systematic. Even if a 
constant radius aperture is used for aperture photometry, 
the error will depend on the sub-pixel shift of the PSF and 
will thus be different for each measured source. Adequate 
sampling is essential, but we believe our tested limit of a > 
1.2 pixels (FWHM > 3.0 pixels) is a reliable benchmark. 



8 SUMMARY 

We have developed an algorithm for performing aperture 
photometry on critically sampled astrophysical images. The 
flux is numerically integrated over the desired aperture, with 
sinc-interpolation used to reconstruct values between pixels. 
As the pixel values can be factored out of the computation- 
ally intensive integral, only the sine functions need to be 
integrated, and these can be pre-computed and applied to 
an arbitrary number of aperture flux measurements. 

The aperture flux is ultimately computed as a weighted 
sum of pixel values, with the sinc-integrals providing the 
weights. The weight values are non-zero outside the nominal 
cut-off radius of the aperture. We find that an aperture 
flux must be measured in a footprint with a border 
extending ~4 pixels beyond the aperture limits to 
ensure accuracy. 

Using Fourier methods, we have shown that the inte- 
gral of the sinc-functions over the aperture is a convolution 
integral, and can be computed quickly as a product in the 
wave-number domain. 

As the method relies on sinc-interpolation based on the 



sampling-theorem, a band-limited PSF is an implicit as- 
sumption. A rule of thumb for this requirement to be satis- 
fied is that a double-Gaussian PSF must be sampled 
such that a > 1.2 pixels, or equivalently FWHM > 3.0 
pixels. No Gaussian-based function is truly band limited, 
but this rule of thumb reduces the photometric error con- 
tributed by the aperture to < 0.001 mag. If the requirement 
is not satisfied, the resulting error depends on the sub-pixel 
position of the PSF and cannot be corrected during calibra- 
tion. The error is also a function of aperture radius, and the 
band-limited requirement must be satisfied to produce an 
accurate growth curve. 

Comparison with SExtractor demonstrated that the 
sine method provides equal performance for large apertures 
(r >5 pixel), and much better performance for smaller aper- 
tures. However, for under-sampled data, the best results 
were obtained with SExtractor using a large aperture. 

The method performs equally well for annular and el- 
liptical apertures, and can be adapted for any aperture. To 
demonstrate the method, a short example has been coded 
in Python and is presented in Appendix [B] 

In closing, we note that the Sloan Digital Sky Survey's 
photometric pipeline used a version of this technique for 
the measurement of aperture fluxes. In that version, the 
Wij components (see equation llOp were pre-computed as 
real-space 2D integrals prior to frame processing, and were 
Lanczos-shifted to be recentred on each source for flux mea- 
surement. As the SDSS photometric catalog is among the 
most widely used data sets in the astronomy community, 
the effectiveness of the algorithm has been very well demon- 
strated in practice. The method is now being used in the 
LSST photometric pipeline under development, with the 
only difference being the use of Fourier methods to more 
efficiently generate the Wij components. Tests on simulated 
LSST images have shown the algorithm to perform effi- 
ciently (i.e., quickly) and produce Poisson-limited accuracy. 



ACKNOWLEDGMENTS 

We wish to thank our LSST colleagues Pat Burchat, Seth 
Digel, Jon Thaler, Phil Marshall, and Rachel Mandelbaum 
for their comments and advice; and we wish to thank Em- 
manuel Bertin for his many valuable suggestions. 

Support for this research was provided by both the 
SDSS and LSST projects. 

Funding for the SDSS and SDSS-II has been provided 
by the Alfred P. Sloan Foundation, the Participating In- 
stitutions, the National Science Foundation, the U.S. De- 
partment of Energy, the National Aeronautics and Space 
Administration, the Japanese Monbukagakusho, the Max 
Planck Society, and the Higher Education Funding Council 
for England. The SDSS Web Site is http://www.sdss.org/ 

The SDSS is managed by the Astrophysical Research 
Consortium for the Participating Institutions. The Partic- 
ipating Institutions are the American Museum of Natu- 
ral History, Astrophysical Institute Potsdam, University of 
Basel, University of Cambridge, Case Western Reserve Uni- 
versity, University of Chicago, Drexel University, Fermilab, 
the Institute for Advanced Study, the Japan Participation 
Group, Johns Hopkins University, the Joint Institute for 
Nuclear Astrophysics, the Kavli Institute for Particle As- 



© 2013 RAS, MNRAS OOO.fnflTI 



10 S.J. Bickerton 



trophysics and Cosmology, the Korean Scientist Group, the 
Chinese Academy of Sciences (LAMOST), Los Alamos Na- 
tional Laboratory, the Max-Planck-Institute for Astronomy 
(MPIA), the Max-Planck-Institute for Astrophysics (MPA), 
New Mexico State University, Ohio State University, Uni- 
versity of Pittsburgh, University of Portsmouth, Princeton 
University, the United States Naval Observatory, and the 
University of Washington. 

LSST project activities are supported in part by the 
National Science Foundation through Governing Coopera- 
tive Agreement 0809409 managed by the Association of Uni- 
versities for Research in Astronomy (AURA), and the De- 
partment of Energy under contract DE-AC02-76-SFO0515 
with the SLAG National Accelerator Laboratory. Additional 
LSST funding comes from private donations, grants to uni- 
versities, and in-kind support from LSSTC Institutional 
Members. 



APPENDIX A: THE FRACTIONAL ALIASED 
POWER IN A SAMPLED GAUSSIAN 

The principal requirement for the sine reconstruction to be 
valid is that the PSF be band-limited. To test this assump- 
tion analytically, we approximate the PSF, p(r), as a circular 
bivariate Gaussian, which has a convenient analytic Hankel 
transform P{q) (circularly symmetric form of the Fourier 
transform) ; 



p[r) — e ' <^ P[q) = a e ^ ' , 



(Al) 



where -i^ denotes a transform pair. 

The aliased power beyond the band-limit for a Gaussian 
(> TT when a is in pixel units) is 



aliased power — 2tv a e 



q dq — Tva e 



(A2) 



where the integrand is the power from equation I All Under 
the assumption that the pixels are square, the convenient 
circular symmetry used here is an approximation, but a very 
good one. The Fourier transform of a real PSF should be 
integrated in two dimensions with limits k:c,y < —T^,kx,y > 
TT, so our use of circular symmetry {q = {k% + ky)^^^ > tt) is 
conservative in that it omits power that is not actually lost 
in Cartesian coordinates (i.e., the corners of the region in 
the wave-number domain where vr < (fc^ -I- ky)^^^). 

Next, we translate the aliased power to a correspond- 
ing error in flux using Rayleigh's Theorerr0 applied to the 
Hankel transform: 



(A3) 



|p(r)| r dr- = / |P(g)| q dq. 



The aliased power equals the integral of the squared 
error |e(r)p, which that power represents in the spatial do- 
main. By taking the square root, we obtain a flux error es- 
timate Q 

^ Rayleigh's Theorem is the continuous form of the (often more 

familiar) discrete summation, Parseval's Theorem. 

* We are interested in the total flux error, and we do not divide by 



error — i'^^ J ^ '^^ 



1/2 



(A4) 



= ( Stt / a'^e " q dq 



1/2 



1/2 -TT^tT''/2 

— n ae 



(A5) 
(A6) 

The total integrated flux is 2Tva^ , and the error shown 
as a fraction is 



total flux 27ri/2cr ' ^ ' 

Equation IA7I can be used as an estimate of the sys- 
tematic flux error introduced by aliased power. Three as- 
sumptions were implicitly made to simplify the estimate: 
(1) circular symmetry can be used as an approximation in 
equation I A21 (2) the power above the Nyquist wave number 
is simply lost rather than being aliased back into the pass 
band, and (3) the flux being integrated is the total flux inte- 
grated to inflnity rather than that within a limited aperture 
radius. All of these assumptions tend to oiierestimate the er- 
ror, and equation I A7I can be regarded as providing an upper 
limit on the systematic flux error. 



APPENDIX B: SAMPLE CODE TO 
DEMONSTRATE SINC APERTURE 
PHOTOMETRY 

The following sample of Python code demonstrates the sine 
method for performing aperture photometry. Python is an 
open-source scripting language, and is readily available on 
most computing platforms (see http://www.python.org). 
Our example uses the numpy and scipy modules (available 
at http : // numpy . scipy . org/^ and http: //www. scipy . org/, 
respectively). 

Due to the nature of Python, the indentation structure 
is essential. The authors will gladly provide a digital copy 
upon request. Note that the code is intended only as an ex- 
ample, and is not written with efflciency in mind. Comput- 
ing performance can be greatly improved by transcribing the 
code to C, C-|— 1-, or some other compiled language. A C-|— I- 
version of the algorithm is included in the LSST code stack. 
For information on downloading the LSST source code see 
www. Isst . org). 

The following command will run the example (assuming 
the script has been named sinc_phot .py), which will plant 
a Gaussian PSF with a — 1.2 centred a,t x,y = 16, 16 in a 
n = 32 X 32 pixel image, and will compute a growth curve 
for it with the sinc-aperture method. Output is printed to 
the screen (i.e., stdout). 

./sinc.phot.py 32 16.0 16.0 1.2 
(i.e., . /sinc_phot .py n x y o") 

The function wi jCoef f icients () is capable of comput- 
ing annular and/or elliptical aperture coefficients, but values 



the aperture area to obtain a per-pixel root- mean-squared (RMS) 
value. The error shown is a root- sum-squared value. 
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for the ellipticity, position angle (theta), and inner radius 
are set to 0.0 in the example. 

#! /usr/bin/ env python 

import sys , numpy , scipy . special as special 

def gaussianlmage (n, center_x, center^y , sigma) : 
image = numpy . zeros C [n, n] , dtype=float) 
for ix in rajigeCn): 

for iy in range (n) : 

X, y = ix - center_x, iy - center.y 
A = 1.0/C2.0*numpy.pi*sigma**2) 

image[iy,ix] = A*numpy.expC-Cx**2 + y**2) / (2 . 0*sigma**2) ) 
return image 

def wij Coefficients (n, x, y , radl , rad2 , theta, ellipticity) : 
wid, xcen, ycen = n, Cn+l)/2, Cn+l)/2 
dx, dy = X - xcen, y - ycen 
if n"/,2: 

dx , dy = dx + 1 , dy + 1 

ftWij = numpy . zeros C [wid, wid] , dtype^complex) 
scale = 1.0 - ellipticity 
for iy in rEingeCwid): 

ky = float (iy - ycen) /wid 

for ix in range (wid): 

kx = float (ix - xcen) /wid 

# rotate ajid rescale 

cosT, sinT = numpy, cos (theta) , numpy. sin(theta) 
kxr, kyr = kx*cosT + ky*sinT, scale*(-kx*sinT + ky*cosT) 
k = numpy . sqrt Ckxr**2 + kyr*+2) 

# compute the airy terms, ajid apply shift theorem 
if k != 0.0: 

airyl = radl* special . j 1(2 . 0*numpy .pi*radl*k) /k 
airy2 = rad2*special . j 1(2 . 0*numpy .pi*rad2*k) /k 
else : 

airyl, airy2 = numpy. pi*radl**2, numpy. pi*rad2**2 
airy = airy2 - airyl 

phase = numpy . exp(-l . 0j*2 . 0*numpy .pi* Cdx*kxr + dy*kyr)) 

ftWij[iy,ix] = phase*scale*airy 

ftWijShift - numpy. fft.fftshift(ftWij) 
wijShift = numpy. fft.ifft2CftWijShift) 
wij = numpy . f ft . fftshift (wijShift) 

return wij . real 

if name ' main ' : 

n = int (sys . argv [1] ) 

X, y, psf_sigma = map(float, sys . argv [2 : ] ) 

radius^inner , theta, ellipticity = 0.0, numpy . pi* (0 . 0) /180 . , 0.0 
psf = gaussianlmage (n, x, y, psf_sigma) 

for radius in numpy. arange(radius_inner + 0.1, 6 . 0*psf _sigma, 0.1): 

wij = wijCoef f icients (n, x, y , radius^inner , radius , theta, ellipticity) 
f lux_measured = (psf *wij ) . sumC) 

f lux.analytic = (1.0 - numpy . exp(-radius**2/ (2 . 0*psf _sigma**2) ) ) 
print radius , f lux_measured, f lux_analytic 
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